Analyses for Dimensions of the language environment in infancy and emerging psychopathology in toddlerhood

Environment

# Libraries
library(tidyverse)
library(lubridate)
library(corrr)
library(corrplot)
library(ggpubr)
library(modelr)
library(effectsize)
library(glmpath)
library(glmnet)
library(MASS)
library(ggcorrplot)
library(broom)
library(ggfortify)
library(relaimpo)
library(olsrr)
library(parameters)
library(BayesFactor)
# for meff function
source('https://osf.io/ws6xc/download')

# Files 
home <- "~/Box/Mooddata_Coordinating/BABIES/Data/final_scored_data/"

eligible_id_file <- "~/Box/lucy_king_files/BABIES/lena_symptoms/inclusion_tracker_sub2.xlsx"
lena_edited_file <- "~/Box/lucy_king_files/BABIES/lena_symptoms/included_recording_dates.csv"
fp_intervention_id <- paste0(home, "lab_caregiving_behavior/free_play_intervention_assignment.csv")
itsea_file <- paste0(home, "survey_18month/itsea_scored_20201001.csv")
demographics_file <- paste0(home, "demographics/demo_6mo_cleaned_final.csv")
lena_file <- paste0(home, "LENA/lena_final_wf_day1_20200208_upto9months.csv")
ibq_file <- paste0(home, "IBQ/ibq_scored_final.csv")
sfp_care_file <- paste0(home, "lab_caregiving_behavior/PCIRS_sfp_wf_complete.csv")
fp_care_file <- paste0(home, "lab_caregiving_behavior/free_play_wf_8min_20200507.csv")
cesd_file <- paste0(home, "CESD/cesd_wf_20201001.csv")

# Functions
source("identify_outliers_histogram.R")

theme_lena <-
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.title = element_text(size = 18), 
    legend.text = element_text(size = 16),
    legend.position = "bottom"
  )

# Display options
options(scipen = 999)

Read in and join data

d0 <-
  readxl::read_xlsx(eligible_id_file, sheet = "data_filtered_T3") %>% 
  dplyr::select(
    ID,
    eligible,
    included
  ) %>% 
  filter(included == 1) %>% 
  left_join(
    read_csv(lena_edited_file) %>% 
      dplyr::select(
        ID, LENA_edited
      ),
    by = "ID"
  ) %>% 
  left_join(
    read_csv(demographics_file) %>% 
      dplyr::select(
        ID,
        baby_dob,
        education,
        education_txt,
        annual_income_txt,
        income_needs,
        male,
        mom_age,
        mom_latinx,
        mom_race,
        age_behav,
        secondlang,
        secondlang_type
      ),
    by = "ID"
  ) %>% 
  left_join(
    read_csv(itsea_file) %>% 
      dplyr::select(
        ID,
        survey_date,
        itsea_symptoms,
        itsea_activity,
        itsea_agress,
        itsea_gad,
        itsea_negemo,
        itsea_depp,
        itsea_intern,
        itsea_extern,
        # itsea symptoms items
        itsea_a2, 
        itsea_a4, 
        itsea_b8, 
        itsea_b10, 
        itsea_b27, 
        itsea_b32,
        itsea_a28, 
        itsea_a30, 
        itsea_a33, 
        itsea_b3, 
        itsea_b9, 
        itsea_b16, 
        itsea_b30, 
        itsea_b34, 
        itsea_b44, 
        itsea_b88, 
        itsea_c2, 
        itsea_e4,
        itsea_a9, 
        itsea_a39, 
        itsea_a40, 
        itsea_b43, 
        itsea_b76, 
        itsea_b81, 
        itsea_b84, 
        itsea_b91, 
        itsea_b92,
        itsea_a3, 
        itsea_a12, 
        itsea_a35, 
        itsea_b37, 
        itsea_b70, 
        itsea_b86, 
        itsea_e5, 
        itsea_e6, 
        itsea_e11,
        itsea_a7, 
        itsea_a21, 
        itsea_a23, 
        itsea_b31, 
        itsea_b45, 
        itsea_b50, 
        itsea_b53, 
        itsea_b59, 
        itsea_b65, 
        itsea_b66, 
        itsea_b74, 
        itsea_b80, 
        itsea_b85,
        itsea_concern_ied
      ), 
    by = "ID"
  ) %>% 
  mutate(
    survey_date =  parse_date_time(survey_date, orders = c("mdy HM")),
    age_18mo = (baby_dob %--% survey_date) / months(1)
  ) %>% 
  left_join(
    read_csv(fp_care_file) %>% 
      dplyr::select(ID, negmood_FP, sens_FP),
    by = "ID"
  ) %>% 
  left_join(
    read_csv(sfp_care_file) %>% 
      dplyr::select(ID, sens_R_M, sens_M),
    by = "ID"
  ) %>% 
  left_join(
    read_csv(ibq_file) %>% 
      dplyr::select(
        ID,
        #IBQ NEG items
        ibq_64, ibq_74, ibq_75, ibq_32, ibq_79, ibq_80,
        ibq_2, ibq_3_r, ibq_4, ibq_21, ibq_52, ibq_53, ibq_62,
        ibq_22, ibq_76, ibq_77, ibq_78, ibq_87, ibq_89,
        ibq_36, ibq_37, ibq_38, ibq_63, ibq_71, ibq_72,
        NEG
      ),
    by = "ID"
  ) %>% 
  left_join(read_csv(lena_file), by = "ID") %>% 
  left_join(read_csv(cesd_file), by = "ID") %>% 
  left_join(read_csv(fp_intervention_id), by = "ID") %>% 
  rename(
    age_lena_d1 = age,
    sens_r_sfp = sens_R_M,
    sens_sfp = sens_M,
    sens_fp = sens_FP,
    negmood_fp = negmood_FP
  ) 

Parsed with column specification:
cols(
  ID = col_double(),
  date_record = col_character(),
  LENA_edited = col_double()
)
Parsed with column specification:
cols(
  .default = col_double(),
  due_date = col_datetime(format = ""),
  baby_dob = col_datetime(format = ""),
  mom_dob = col_datetime(format = ""),
  baby_race = col_character(),
  baby_race_describe = col_character(),
  mom_race = col_character(),
  momrace_describe = col_character(),
  annual_income_txt = col_character(),
  education_txt = col_character(),
  employment_status_txt = col_character(),
  marital_status_txt = col_character(),
  partner_educ_txt = col_character(),
  partner_employ_txt = col_character(),
  questionnaire_only_date = col_date(format = ""),
  educ_describe = col_logical(),
  employment_explain = col_character(),
  mom_pob = col_character(),
  mom_native_lang = col_character(),
  primarylang = col_character(),
  secondlang_type = col_character()
  # ... with 24 more columns
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  survey_date = col_character(),
  baby_dob = col_datetime(format = ""),
  activity_concern = col_character(),
  agress_concern = col_character(),
  depp_concern = col_character(),
  gad_concern = col_character(),
  negemo_concern = col_character(),
  play_concern = col_character(),
  empathy_concern = col_character(),
  social_concern = col_character(),
  itsea_concern = col_character(),
  itsea_concern_ied = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  ID = col_double(),
  sens_FP = col_double(),
  intrus_FP = col_double(),
  posreg_FP = col_double(),
  stim_FP = col_double(),
  negreg_FP = col_double(),
  negmood_FP = col_double(),
  posmood_FP = col_double(),
  detach_FP = col_double()
)
Parsed with column specification:
cols(
  .default = col_double()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ibq_timestamp = col_datetime(format = "")
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  date_record = col_date(format = ""),
  day_type = col_character(),
  lena_recordno = col_character(),
  lena_probyes = col_character(),
  other_caregiver = col_character(),
  lena_notes = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  ID = col_double(),
  cesd_t1 = col_double(),
  cesd_t2 = col_double(),
  cesd_t3 = col_double(),
  cesd_t4 = col_double(),
  cesd_t1_concern = col_character(),
  cesd_t2_concern = col_character(),
  cesd_t3_concern = col_character(),
  cesd_t4_concern = col_character()
)
Parsed with column specification:
cols(
  ID = col_double(),
  group = col_character()
)

Retention / Missing data

Dyads who completed the 6-month assessment (either lab or home observation), were eligible for the 18-month assessment (i.e., were age 18 months by the time of analysis), and completed the 18-month assessment are included in the current analyses.

#missing SFP
d0 %>% 
  count(is.na(sens_sfp))

#missing FP
d0 %>% 
  count(is.na(sens_fp))

#missing IBQ
d0 %>% 
  count(is.na(NEG))

d0 %>% 
  dplyr::select(
    awc_prop, 
    ctc_prop,
    awc_hour_max,
    ctc_hour_max,
    sens_r_sfp,
    sens_fp,
    NEG,
  ) %>% 
  pair_n()

             awc_prop ctc_prop awc_hour_max ctc_hour_max sens_r_sfp sens_fp NEG
awc_prop          100      100          100          100         94      96  96
ctc_prop          100      100          100          100         94      96  96
awc_hour_max      100      100          100          100         94      96  96
ctc_hour_max      100      100          100          100         94      96  96
sens_r_sfp         94       94           94           94         94      93  90
sens_fp            96       96           96           96         93      96  92
NEG                96       96           96           96         90      92  96
attr(,"class")
[1] "n_mat"  "matrix"

Descriptive statistics

IBQ, ITSEA reliability

d0 %>% 
  dplyr::select(ibq_64:ibq_72) %>% 
  psych::alpha(check.keys = TRUE)

Some items were negatively correlated with total scale and were automatically reversed.
 This is indicated by a negative sign for the variable name.

Reliability analysis   
Call: psych::alpha(x = ., check.keys = TRUE)

 

 lower alpha upper     95% confidence boundaries
0.84 0.87 0.91 

 Reliability if an item is dropped:

 Item statistics 

Non missing response frequency for each item
           1    2    3    4    5    6    7 miss
ibq_64  0.07 0.20 0.16 0.12 0.12 0.20 0.13 0.10
ibq_74  0.04 0.12 0.13 0.16 0.23 0.24 0.07 0.06
ibq_75  0.29 0.35 0.11 0.06 0.10 0.08 0.01 0.11
ibq_32  0.44 0.45 0.06 0.03 0.00 0.01 0.00 0.07
ibq_79  0.08 0.21 0.13 0.28 0.09 0.17 0.04 0.08
ibq_80  0.03 0.24 0.20 0.23 0.11 0.18 0.01 0.07
ibq_2   0.05 0.17 0.20 0.27 0.12 0.16 0.03 0.05
ibq_3_r 0.01 0.18 0.19 0.32 0.09 0.18 0.03 0.05
ibq_4   0.00 0.22 0.17 0.22 0.16 0.18 0.06 0.04
ibq_21  0.02 0.23 0.25 0.19 0.16 0.12 0.03 0.04
ibq_52  0.02 0.09 0.22 0.18 0.17 0.26 0.05 0.05
ibq_53  0.34 0.27 0.14 0.06 0.12 0.04 0.03 0.06
ibq_62  0.04 0.22 0.14 0.21 0.11 0.20 0.07 0.10
ibq_22  0.07 0.25 0.11 0.17 0.17 0.18 0.06 0.11
ibq_76  0.24 0.30 0.17 0.08 0.05 0.10 0.06 0.07
ibq_77  0.47 0.26 0.06 0.06 0.05 0.08 0.03 0.12
ibq_78  0.51 0.30 0.08 0.06 0.02 0.02 0.01 0.11
ibq_87  0.61 0.26 0.09 0.00 0.04 0.00 0.00 0.11
ibq_89  0.47 0.31 0.07 0.03 0.01 0.08 0.03 0.28
ibq_36  0.00 0.07 0.23 0.15 0.18 0.33 0.04 0.05
ibq_37  0.01 0.40 0.29 0.16 0.06 0.04 0.03 0.05
ibq_38  0.02 0.04 0.12 0.17 0.18 0.45 0.02 0.04
ibq_63  0.01 0.09 0.11 0.19 0.22 0.37 0.01 0.05
ibq_71  0.00 0.02 0.01 0.11 0.16 0.42 0.28 0.07
ibq_72  0.55 0.37 0.05 0.00 0.01 0.01 0.00 0.06

d0 %>% 
  dplyr::select(itsea_a2:itsea_b85, -itsea_c2) %>% 
  psych::alpha()

Item = itsea_e4 had no variance and was deleted but still is counted in the scoreSome items were negatively correlated with the total scale and probably 
should be reversed.  
To do this, run the function again with the 'check.keys=TRUE' option

Some items ( itsea_b84 itsea_a35 itsea_b37 itsea_e5 ) were negatively correlated with the total scale and 
probably should be reversed.  
To do this, run the function again with the 'check.keys=TRUE' option
Reliability analysis   
Call: psych::alpha(x = .)

 

 lower alpha upper     95% confidence boundaries
0.79 0.83 0.88 

 Reliability if an item is dropped:

 Item statistics 

Non missing response frequency for each item
             0    1    2 miss
itsea_a2  0.70 0.29 0.01 0.00
itsea_a4  0.46 0.41 0.13 0.00
itsea_b8  0.31 0.51 0.18 0.00
itsea_b10 0.19 0.56 0.25 0.00
itsea_b27 0.55 0.39 0.06 0.00
itsea_b32 0.82 0.14 0.04 0.01
itsea_a28 0.80 0.19 0.01 0.00
itsea_a30 0.72 0.28 0.00 0.00
itsea_a33 0.90 0.09 0.01 0.00
itsea_b3  0.53 0.45 0.02 0.00
itsea_b9  0.64 0.32 0.04 0.00
itsea_b16 0.67 0.31 0.02 0.00
itsea_b30 0.47 0.51 0.02 0.00
itsea_b34 0.92 0.07 0.01 0.00
itsea_b44 0.19 0.63 0.18 0.00
itsea_b88 0.41 0.54 0.05 0.00
itsea_a9  0.97 0.03 0.00 0.00
itsea_a39 0.94 0.06 0.00 0.00
itsea_a40 0.98 0.02 0.00 0.00
itsea_b43 0.78 0.20 0.02 0.01
itsea_b76 0.96 0.04 0.00 0.02
itsea_b81 0.97 0.03 0.00 0.00
itsea_b84 0.98 0.01 0.01 0.00
itsea_b91 0.99 0.01 0.00 0.00
itsea_b92 0.98 0.02 0.00 0.00
itsea_a3  0.89 0.10 0.01 0.00
itsea_a12 0.83 0.15 0.02 0.01
itsea_a35 0.80 0.17 0.03 0.00
itsea_b37 0.98 0.02 0.00 0.00
itsea_b70 0.78 0.19 0.03 0.01
itsea_b86 0.72 0.25 0.03 0.00
itsea_e5  0.77 0.21 0.02 0.00
itsea_e6  0.81 0.18 0.01 0.00
itsea_e11 0.92 0.07 0.01 0.00
itsea_a7  0.72 0.23 0.05 0.00
itsea_a21 0.58 0.40 0.02 0.01
itsea_a23 0.60 0.35 0.05 0.00
itsea_b31 0.10 0.77 0.13 0.00
itsea_b45 0.77 0.22 0.01 0.00
itsea_b50 0.69 0.27 0.04 0.00
itsea_b53 0.72 0.28 0.00 0.00
itsea_b59 0.42 0.51 0.07 0.00
itsea_b65 0.11 0.73 0.16 0.00
itsea_b66 0.82 0.16 0.02 0.01
itsea_b74 0.77 0.23 0.00 0.02
itsea_b80 0.78 0.17 0.05 0.01
itsea_b85 0.55 0.44 0.01 0.00

d0 %>% 
  summarize_at(
    vars(
      age_18mo,
      age_lena_d1,
      NEG,
      duration_total,
      duration_analyzed,
      ct_first_time,
      ct_last_time,
      cvc_hour_max,
      ctc_hour_max,
      awc_hour_max,
      ctc_prop,
      awc_prop,
      sens_r_sfp,
      sens_fp,
      negmood_fp,
      itsea_symptoms,
      mom_age,
      income_needs,
      percent_mother,
      cesd_t3,
      cesd_t4
    ), 
    funs(mean, sd, min, max), na.rm = TRUE
  )

`funs()` is deprecated as of dplyr 0.8.0.
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

d0 %>% 
  count(male) %>% 
  mutate(
    per = n / sum(n)
  )

d0 %>% 
  count(itsea_concern_ied) %>% 
  mutate(
    per = n / sum(n)
  )

d0 %>% 
  count(income_needs < 1) %>% 
  mutate(
    per = n / sum(n)
  )

d0 %>% 
  count(mom_latinx) %>% 
  mutate(
    per = n / sum(n)
  )

d0 %>% 
  count(mom_race) %>% 
  mutate(
    per = n / sum(n)
  ) %>% 
  arrange(desc(n))

d0 %>% 
  count(education_txt) %>% 
  mutate(
    per = n / sum(n)
  ) %>% 
  arrange(education_txt)

d0 %>% 
  count(day_type) %>% 
  mutate(per = n / sum(n))

d0 %>% 
  count(cesd_t3 >= 16) %>% 
  mutate(per = n / sum(n))

d0 %>% 
  count(cesd_t4 >= 16) %>% 
  mutate(per = n / sum(n))

NA

d0 %>% 
  count(lena_mornstart, lena_recordfull)

d0 %>% 
  filter(lena_mornstart == 0) %>% 
  dplyr::select(
    ID, 
    lena_recordno, 
    ct_first_time, 
    ct_last_time, 
    duration_total, 
    duration_analyzed, 
    lena_notes
  )

Distributions

d0 %>% 
  identify_outliers_hist(x = cvc_hour_max)

d0 %>% 
  identify_outliers_hist(x = ctc_hour_max)

d0 %>% 
  identify_outliers_hist(x = ctc_prop)

d0 %>% 
  identify_outliers_hist(x = awc_hour_max)

d0 %>% 
  identify_outliers_hist(x = awc_prop)

d0 %>% 
  identify_outliers_hist(x = sens_r_sfp)

d0 %>% 
  identify_outliers_hist(x = sens_fp)

d0 %>% 
  identify_outliers_hist(x = itsea_symptoms)

d0 %>% 
  identify_outliers_hist(x = cesd_t3)

d0 %>% 
  identify_outliers_hist(x = cesd_t4)

d0 %>% 
  ggplot(aes(percent_mother)) +
  geom_histogram() +
  theme_lena +
  labs(
    x = "% of the LENA recording day infant with mother"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/percent_mom_hist.png",
  width = 7,
  height = 5
)

Correlations

corr_lena <- 
  d0 %>%
  dplyr::select(
    `AW consistency` = awc_prop,
    `CT consistency` = ctc_prop,
    `AW quantity` = awc_hour_max,
    `CT quantity` = ctc_hour_max,
    `Infant vocalizations` = cvc_hour_max,
    `Sensitivity to distress` = sens_r_sfp,
    `Sensitivity during play` = sens_fp,
  ) %>% 
  correlate(use = "pairwise.complete.obs", method = "pearson") %>% 
  shave(upper = TRUE) %>% 
  fashion()  

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

corr_lena_plot <-
  corr_lena %>% 
  gather(variable, value, -rowname) %>% 
  mutate(
    variable = str_replace_all(variable, "\\.", " "),
    rowname = as.character(rowname),
    value_chr = as.character(value),
    value_chr = case_when(
      rowname == "CT consistency" & variable == "AW quantity" ~ ".23",
      rowname == "AW quantity" & variable == "CT consistency" ~ NA_character_,
      rowname == "Sensitivity to distress" & variable == "Sensitivity during play" ~ ".10",
      rowname == "Sensitivity during play" & variable == "Sensitivity to distress" ~ NA_character_,
      TRUE ~ value_chr
    ),
    value_num = as.numeric(value_chr)
  ) 

corr_lena_plot %>% 
  ggplot(aes(x = rowname, y = variable)) +
  geom_tile(aes(fill = abs(value_num))) +
  geom_text(
    aes(label = value_chr), 
    size = 4
  ) +
  scale_fill_gradient2(
    low = "blue",
    high = "red",
    mid = "white",
    na.value = "white",
    midpoint = 0,
    limit = c(-1, 1),
    space = "Lab",
    name = "Pearson correlation\ncoefficient"
  ) +
  theme_void() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1,
      size = 12
    ),
    axis.text.y = element_text(
      size = 12,
      hjust = 1.1
    )
  ) 

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/lena_corr_plot.eps", 
  dpi = 1000,
  height = 5, 
  width = 7
)

cor.test(d0$cesd_t3, d0$awc_hour_max)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$awc_hour_max
t = 0.70855, df = 98, p-value = 0.4803
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1268048  0.2641061
sample estimates:
       cor 
0.07139142 

cor.test(d0$cesd_t3, d0$awc_prop)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$awc_prop
t = 0.20234, df = 98, p-value = 0.8401
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1766922  0.2159863
sample estimates:
       cor 
0.02043509 

cor.test(d0$cesd_t3, d0$ctc_hour_max)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$ctc_hour_max
t = 0.30729, df = 98, p-value = 0.7593
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1664066  0.2260661
sample estimates:
       cor 
0.03102564 

cor.test(d0$cesd_t3, d0$ctc_prop)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$ctc_prop
t = -0.040919, df = 98, p-value = 0.9674
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2003889  0.1924409
sample estimates:
        cor 
-0.00413343 

cor.test(d0$cesd_t3, d0$cvc_hour_max)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$cvc_hour_max
t = 0.046727, df = 98, p-value = 0.9628
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1918759  0.2009519
sample estimates:
       cor 
0.00472013 

cor.test(d0$cesd_t3, d0$sens_r_sfp)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$sens_r_sfp
t = 0.75743, df = 92, p-value = 0.4507
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1259026  0.2769225
sample estimates:
       cor 
0.07872259 

cor.test(d0$cesd_t3, d0$sens_fp)

    Pearson's product-moment correlation

data:  d0$cesd_t3 and d0$sens_fp
t = -0.18924, df = 94, p-value = 0.8503
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2191435  0.1816818
sample estimates:
        cor 
-0.01951496 

cor.test(d0$NEG, d0$awc_hour_max)

    Pearson's product-moment correlation

data:  d0$NEG and d0$awc_hour_max
t = 0.81245, df = 94, p-value = 0.4186
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1189726  0.2793149
sample estimates:
       cor 
0.08350513 

cor.test(d0$NEG, d0$awc_prop)

    Pearson's product-moment correlation

data:  d0$NEG and d0$awc_prop
t = -0.52757, df = 94, p-value = 0.599
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2520746  0.1477609
sample estimates:
        cor 
-0.05433463 

cor.test(d0$NEG, d0$ctc_hour_max)

    Pearson's product-moment correlation

data:  d0$NEG and d0$ctc_hour_max
t = 1.4766, df = 94, p-value = 0.1431
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.05147543  0.34076399
sample estimates:
      cor 
0.1505644 

cor.test(d0$NEG, d0$ctc_prop)

    Pearson's product-moment correlation

data:  d0$NEG and d0$ctc_prop
t = -0.3159, df = 94, p-value = 0.7528
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2315400  0.1690237
sample estimates:
        cor 
-0.03256579 

cor.test(d0$NEG, d0$cvc_hour_max)

    Pearson's product-moment correlation

data:  d0$NEG and d0$cvc_hour_max
t = 0.24743, df = 94, p-value = 0.8051
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1758738  0.2248476
sample estimates:
       cor 
0.02551173 

cor.test(d0$NEG, d0$sens_r_sfp)

    Pearson's product-moment correlation

data:  d0$NEG and d0$sens_r_sfp
t = -0.2727, df = 88, p-value = 0.7857
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2347363  0.1791112
sample estimates:
        cor 
-0.02905776 

cor.test(d0$NEG, d0$sens_fp)

    Pearson's product-moment correlation

data:  d0$NEG and d0$sens_fp
t = -1.1946, df = 90, p-value = 0.2354
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3215266  0.0819784
sample estimates:
       cor 
-0.1249368 

cor.test(d0$awc_hour_max, d0$sens_fp)
cor.test(d0$awc_hour_max, d0$sens_r_sfp)

cor.test(d0$ctc_hour_max, d0$sens_fp)
cor.test(d0$ctc_hour_max, d0$sens_r_sfp)

cor.test(d0$awc_prop, d0$sens_fp)
cor.test(d0$awc_prop, d0$sens_r_sfp)

cor.test(d0$ctc_prop, d0$sens_fp)
cor.test(d0$ctc_prop, d0$sens_r_sfp)

cor.test(d0$cvc_hour_max, d0$sens_fp)
cor.test(d0$cvc_hour_max, d0$sens_r_sfp)

cor.test(d0$sens_fp, d0$sens_r_sfp)

t-tests by bilingual status

t.test(d0$awc_hour_max ~ d0$secondlang)

    Welch Two Sample t-test

data:  d0$awc_hour_max by d0$secondlang
t = 0.28166, df = 94.674, p-value = 0.7788
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -503.7566  670.3237
sample estimates:
mean in group 0 mean in group 1 
       3549.149        3465.865 

t.test(d0$ctc_hour_max ~ d0$secondlang)

    Welch Two Sample t-test

data:  d0$ctc_hour_max by d0$secondlang
t = 0.47747, df = 91.869, p-value = 0.6342
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.311072 15.204689
sample estimates:
mean in group 0 mean in group 1 
       75.44681        72.50000 

t.test(d0$ctc_prop ~ d0$secondlang)

    Welch Two Sample t-test

data:  d0$ctc_prop by d0$secondlang
t = 0.87779, df = 96.693, p-value = 0.3822
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.02600516  0.06724637
sample estimates:
mean in group 0 mean in group 1 
      0.5561746       0.5355540 

t.test(d0$awc_prop ~ d0$secondlang)

    Welch Two Sample t-test

data:  d0$awc_prop by d0$secondlang
t = 0.52872, df = 96.086, p-value = 0.5982
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.03152325  0.05441354
sample estimates:
mean in group 0 mean in group 1 
      0.7422157       0.7307705 

t.test(d0$itsea_symptoms ~ d0$secondlang)

    Welch Two Sample t-test

data:  d0$itsea_symptoms by d0$secondlang
t = -0.34477, df = 96.903, p-value = 0.731
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2353910  0.1657158
sample estimates:
mean in group 0 mean in group 1 
       1.053025        1.087863 

t-tests by day type

t.test(d0$awc_hour_max ~ d0$day_type)

    Welch Two Sample t-test

data:  d0$awc_hour_max by d0$day_type
t = -1.4866, df = 66.225, p-value = 0.1419
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1072.0120   156.9219
sample estimates:
mean in group weekday mean in group weekend 
             3348.369              3805.914 

t.test(d0$ctc_hour_max ~ d0$day_type)

    Welch Two Sample t-test

data:  d0$ctc_hour_max by d0$day_type
t = 1.0076, df = 69.607, p-value = 0.3171
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.267328 19.062933
sample estimates:
mean in group weekday mean in group weekend 
             76.36923              69.97143 

t.test(d0$ctc_prop ~ d0$day_type)

    Welch Two Sample t-test

data:  d0$ctc_prop by d0$day_type
t = 0.74512, df = 65.551, p-value = 0.4589
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.03151406  0.06903400
sample estimates:
mean in group weekday mean in group weekend 
            0.5534388             0.5346789 

t.test(d0$awc_prop ~ d0$day_type)

    Welch Two Sample t-test

data:  d0$awc_prop by d0$day_type
t = -1.4395, df = 68.83, p-value = 0.1545
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.07686127  0.01243104
sample estimates:
mean in group weekday mean in group weekend 
            0.7255763             0.7577915 

Histograms

d0 %>% 
  dplyr::select(
    ctc_hour_max,
    awc_hour_max,
    ctc_prop,
    awc_prop,
    sens_fp,
    sens_r_sfp
  ) %>% 
  gather(measure, value, ctc_hour_max:sens_r_sfp) %>% 
  mutate(
    measure = factor(
      measure,
      levels = c(
        "awc_prop",
        "ctc_prop",
        "awc_hour_max",
        "ctc_hour_max",
        "sens_r_sfp",
        "sens_fp"
      ),
      labels = c(
        "AW consistency", 
        "CT consistency", 
        "AW quantity", 
        "CT quantity",
        "Sensitivity to distress",
        "Sensivity during play"
        )
    )
  ) %>% 
  ggplot(aes(value)) +
  geom_histogram() +
  facet_wrap(.~measure, scales = "free", ncol = 2, nrow = 3) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 18),
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.position = "none"
  ) +
  labs(
    x = "Caregiving input value"
  )


ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/care_hist_ggpanels.png",
  width = 7.5,
  height = 8
)

d0 %>% 
  ggplot(aes(itsea_symptoms)) +
  geom_histogram() +
  scale_x_continuous(breaks = seq.int(0, 3, .25)) +
  theme_lena +
  labs(
    x = "Toddler symptoms of psychopathology"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/itsea_hist.png"
)

Saving 6.33 x 3.91 in image

d0 %>% 
  dplyr::select(
    Activity = itsea_activity,
    Aggression = itsea_agress,
    Depression = itsea_depp,
    Anxiety = itsea_gad,
    `Negative emotion` = itsea_negemo
  ) %>% 
  gather(itsea_subscale, score, Activity:`Negative emotion`) %>% 
  ggplot(aes(score)) +
  geom_histogram(bins = 12) +
  theme_lena +
  theme(
    strip.text = element_text(size = 16)
  ) +
  facet_wrap(.~itsea_subscale, scales = "free") +
  labs(
    x = "Subscale score"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/itsea_subscale_hist.png",
  width = 10,
  height = 7
)

d0 %>% 
  ggplot(aes(cesd_t3)) +
  geom_histogram() +
  scale_x_continuous(breaks = seq.int(0, 50, 5)) +
  theme_lena +
  labs(
    x = "Maternal depressive symptoms\n at infant age 6 months"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cesd6_hist.png"
)

Saving 6.33 x 3.91 in image

d0 %>% 
  ggplot(aes(cesd_t4)) +
  geom_histogram() +
  scale_x_continuous(breaks = seq.int(0, 50, 5)) +
  theme_lena +
  labs(
    x = "Maternal depressive symptoms\n at infant age 18 months"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cesd18_hist.png"
)

Saving 6.33 x 3.91 in image

Results

Associations of LENA measures with ITSEA symptoms

Effect code categorical variables

d0 <-
  d0 %>% 
  mutate(
    white = as.factor(
      if_else(
        mom_race == "White", 
        1, 0
      )
    ),
    male = as.factor(male),
    mom_latinx = as.factor(mom_latinx)
  )

contrasts(d0$male) <- c(-.5, .5)
contrasts(d0$white) <- c(-.5, .5)
contrasts(d0$mom_latinx) <- c(-.5, .5)

Formal model fitting

d0_lm <-
  d0 %>% 
  filter(!is.na(NEG))

mod_0 <-
  lm(
    scale(itsea_symptoms) ~
      1,
    data = d0_lm
  )

mod_cesdt3 <-
  lm(
    scale(itsea_symptoms) ~
      scale(cesd_t3),
    data = d0_lm
  )

anova(mod_0, mod_cesdt3)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ 1
Model 2: scale(itsea_symptoms) ~ scale(cesd_t3)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     95 95.000                           
2     94 93.846  1    1.1538 1.1557 0.2851

mod_cesdt4 <-
  lm(
    scale(itsea_symptoms) ~
      scale(cesd_t4),
    data = d0_lm
  )

anova(mod_0, mod_cesdt4)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ 1
Model 2: scale(itsea_symptoms) ~ scale(cesd_t4)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)  
1     95 95.000                             
2     94 91.346  1    3.6537 3.7598 0.0555 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mod_neg <-
  lm(
    scale(itsea_symptoms) ~
      scale(cesd_t4) +
      scale(NEG),
    data = d0_lm
  )

anova(mod_cesdt4, mod_neg)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4)
Model 2: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1     94 91.346                                 
2     93 81.119  1    10.227 11.725 0.000919 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

model_parameters(mod_neg)

Parameter   | Coefficient |   SE |        95% CI |        t | df |      p
-------------------------------------------------------------------------
(Intercept) |    9.39e-17 | 0.10 | [-0.19, 0.19] | 9.85e-16 | 93 | > .999
cesd_t4     |        0.18 | 0.10 | [-0.01, 0.37] |     1.84 | 93 | 0.069 
NEG         |        0.33 | 0.10 | [ 0.14, 0.52] |     3.42 | 93 | < .001

calc.relimp(mod_neg)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

2 Regressors: 
scale(cesd_t4) scale(NEG) 
Proportion of variance explained by model: 14.61%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                      lmg
scale(cesd_t4) 0.03476271
scale(NEG)     0.11135116

Average coefficients for different model sizes: 

                      1X       2Xs
scale(cesd_t4) 0.1961112 0.1765674
scale(NEG)     0.3391873 0.3286887

mod_Bage <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(age_lena_d1),
    data = d0_lm
  )

anova(mod_neg, mod_Bage)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(age_lena_d1)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 79.246  1    1.8733 2.1748 0.1437

mod_Mage <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(mom_age),
    data = d0_lm
  )

anova(mod_neg, mod_Mage)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(mom_age)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 78.782  1    2.3375 2.7298 0.1019

mod_race <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      white, 
    data = d0_lm
  )

anova(mod_neg, mod_race)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + white
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 80.305  1   0.81468 0.9333 0.3365

mod_ethnicity <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      mom_latinx, 
    data = d0_lm
  )

anova(mod_neg, mod_ethnicity)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + mom_latinx
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 79.287  1    1.8318 2.1255 0.1483

mod_education <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(education),
    data = d0_lm
  )

anova(mod_neg, mod_education)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(education)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 79.548  1    1.5715 1.8175 0.1809

mod_income_needs <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(income_needs), 
    data = d0_lm
  )

anova(mod_neg, mod_income_needs)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(income_needs)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 80.995  1   0.12377 0.1406 0.7086

mod_sex <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      male, 
    data = d0_lm
  )

anova(mod_neg, mod_sex)

Analysis of Variance Table

Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + male
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     93 81.119                           
2     92 79.115  1     2.004 2.3304 0.1303

Consistency

AW

mod_awprop <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(awc_prop),
    data = d0_lm
  )

model_parameters(mod_awprop)

Parameter   | Coefficient |   SE |         95% CI |        t | df |      p
--------------------------------------------------------------------------
(Intercept) |    1.22e-16 | 0.09 | [-0.19,  0.19] | 1.31e-15 | 92 | > .999
NEG         |        0.32 | 0.09 | [ 0.13,  0.50] |     3.36 | 92 | 0.001 
cesd_t4     |        0.20 | 0.09 | [ 0.01,  0.38] |     2.09 | 92 | 0.040 
awc_prop    |       -0.22 | 0.09 | [-0.41, -0.03] |    -2.32 | 92 | 0.023 

calc.relimp(mod_awprop)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

3 Regressors: 
scale(NEG) scale(cesd_t4) scale(awc_prop) 
Proportion of variance explained by model: 19.33%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                       lmg
scale(NEG)      0.10717209
scale(cesd_t4)  0.03851538
scale(awc_prop) 0.04762374

Average coefficients for different model sizes: 

                        1X        2Xs        3Xs
scale(NEG)       0.3391873  0.3284913  0.3156267
scale(cesd_t4)   0.1961112  0.1967953  0.1965872
scale(awc_prop) -0.2183248 -0.2189630 -0.2184895

CT

mod_ctprop <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(ctc_prop),
    data = d0_lm
  )

model_parameters(mod_ctprop)

Parameter   | Coefficient |   SE |         95% CI |        t | df |      p
--------------------------------------------------------------------------
(Intercept) |    1.63e-17 | 0.09 | [-0.18,  0.18] | 1.76e-16 | 92 | > .999
NEG         |        0.32 | 0.09 | [ 0.13,  0.50] |     3.43 | 92 | < .001
cesd_t4     |        0.20 | 0.09 | [ 0.02,  0.39] |     2.16 | 92 | 0.034 
ctc_prop    |       -0.25 | 0.09 | [-0.43, -0.06] |    -2.63 | 92 | 0.010 

calc.relimp(mod_ctprop)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

3 Regressors: 
scale(NEG) scale(cesd_t4) scale(ctc_prop) 
Proportion of variance explained by model: 20.6%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                       lmg
scale(NEG)      0.10841060
scale(cesd_t4)  0.03952569
scale(ctc_prop) 0.05808721

Average coefficients for different model sizes: 

                        1X        2Xs        3Xs
scale(NEG)       0.3391873  0.3302660  0.3191677
scale(cesd_t4)   0.1961112  0.1993232  0.2018538
scale(ctc_prop) -0.2363189 -0.2420650 -0.2461931

Quantity

AW

mod_awmax<-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(awc_hour_max),
    data = d0_lm
  )

model_parameters(mod_awmax)

Parameter    | Coefficient |   SE |        95% CI |        t | df |      p
--------------------------------------------------------------------------
(Intercept)  |    8.29e-17 | 0.10 | [-0.19, 0.19] | 8.67e-16 | 92 | > .999
NEG          |        0.33 | 0.10 | [ 0.14, 0.53] |     3.45 | 92 | < .001
cesd_t4      |        0.18 | 0.10 | [-0.01, 0.37] |     1.86 | 92 | 0.066 
awc_hour_max |       -0.06 | 0.10 | [-0.25, 0.13] |    -0.62 | 92 | 0.537 

calc.relimp(mod_awmax)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

3 Regressors: 
scale(NEG) scale(cesd_t4) scale(awc_hour_max) 
Proportion of variance explained by model: 14.97%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                            lmg
scale(NEG)          0.112514996
scale(cesd_t4)      0.035192701
scale(awc_hour_max) 0.001947467

Average coefficients for different model sizes: 

                             1X         2Xs         3Xs
scale(NEG)           0.33918734  0.33605409  0.33350200
scale(cesd_t4)       0.19611119  0.18725146  0.17959080
scale(awc_hour_max) -0.02200461 -0.04182114 -0.05979412

CT

mod_ctmax <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(ctc_hour_max),
    data = d0_lm
  )

model_parameters(mod_ctmax)

Parameter    | Coefficient |   SE |         95% CI |        t | df |      p
---------------------------------------------------------------------------
(Intercept)  |    5.46e-17 | 0.09 | [-0.19,  0.19] | 5.86e-16 | 92 | > .999
NEG          |        0.36 | 0.09 | [ 0.17,  0.55] |     3.79 | 92 | < .001
cesd_t4      |        0.21 | 0.10 | [ 0.02,  0.40] |     2.23 | 92 | 0.028 
ctc_hour_max |       -0.22 | 0.10 | [-0.41, -0.03] |    -2.28 | 92 | 0.025 

calc.relimp(mod_ctmax)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

3 Regressors: 
scale(NEG) scale(cesd_t4) scale(ctc_hour_max) 
Proportion of variance explained by model: 19.19%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                           lmg
scale(NEG)          0.12030047
scale(cesd_t4)      0.04075023
scale(ctc_hour_max) 0.03087917

Average coefficients for different model sizes: 

                            1X        2Xs        3Xs
scale(NEG)           0.3391873  0.3478042  0.3596103
scale(cesd_t4)       0.1961112  0.2007037  0.2123774
scale(ctc_hour_max) -0.1289444 -0.1758479 -0.2195135

Infant vocalizations

mod_cvmax <-
  lm(
    scale(itsea_symptoms) ~
      scale(NEG) +
      scale(cesd_t4) +
      scale(cvc_hour_max),
    data = d0_lm
  )

model_parameters(mod_cvmax)

Parameter    | Coefficient |   SE |         95% CI |        t | df |      p
---------------------------------------------------------------------------
(Intercept)  |    5.25e-17 | 0.09 | [-0.18,  0.18] | 5.67e-16 | 92 | > .999
NEG          |        0.33 | 0.09 | [ 0.15,  0.52] |     3.56 | 92 | < .001
cesd_t4      |        0.22 | 0.09 | [ 0.03,  0.40] |     2.28 | 92 | 0.025 
cvc_hour_max |       -0.24 | 0.09 | [-0.43, -0.05] |    -2.54 | 92 | 0.013 

calc.relimp(mod_cvmax)

Response variable: scale(itsea_symptoms) 
Total response variance: 1 
Analysis based on 96 observations 

3 Regressors: 
scale(NEG) scale(cesd_t4) scale(cvc_hour_max) 
Proportion of variance explained by model: 20.21%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                           lmg
scale(NEG)          0.11276570
scale(cesd_t4)      0.04185695
scale(cvc_hour_max) 0.04746378

Average coefficients for different model sizes: 

                            1X        2Xs        3Xs
scale(NEG)           0.3391873  0.3365568  0.3325144
scale(cesd_t4)       0.1961112  0.2052806  0.2151002
scale(cvc_hour_max) -0.1965144 -0.2198211 -0.2397697

Bayes Factors

model_d <-
  d0_lm %>% 
  dplyr::select(
    awc_prop,
    ctc_prop,
    awc_hour_max,
    ctc_hour_max,
    cvc_hour_max,
    itsea_symptoms,
    NEG,
    cesd_t4
  ) %>% 
  as.data.frame()

# each measure is compared to a covariate only model
lmBF_covonly = lmBF(itsea_symptoms ~ NEG + cesd_t4, data = model_d)

lmBF_awcprop = lmBF(itsea_symptoms ~ awc_prop + cesd_t4 + NEG, data = model_d)
lmBF_ctcprop = lmBF(itsea_symptoms ~ ctc_prop + cesd_t4 + NEG, data = model_d)

lmBF_awcmax = lmBF(itsea_symptoms ~ awc_hour_max + cesd_t4 + NEG, data = model_d)
lmBF_ctcmax = lmBF(itsea_symptoms ~ ctc_hour_max + cesd_t4 + NEG, data = model_d)

mBF_cvcmax = lmBF(itsea_symptoms ~ cvc_hour_max + cesd_t4 + NEG, data = model_d)

# BF for AW consistency
lmBF_awcprop / lmBF_covonly

Bayes factor analysis
--------------
[1] awc_prop + cesd_t4 + NEG : 3.008407 ±0%

Against denominator:
  itsea_symptoms ~ NEG + cesd_t4 
---
Bayes factor type: BFlinearModel, JZS

# BF for CT consistency
lmBF_ctcprop / lmBF_covonly

Bayes factor analysis
--------------
[1] ctc_prop + cesd_t4 + NEG : 5.794294 ±0%

Against denominator:
  itsea_symptoms ~ NEG + cesd_t4 
---
Bayes factor type: BFlinearModel, JZS

# BF for AW quantity
lmBF_awcmax / lmBF_covonly

Bayes factor analysis
--------------
[1] awc_hour_max + cesd_t4 + NEG : 0.3489947 ±0%

Against denominator:
  itsea_symptoms ~ NEG + cesd_t4 
---
Bayes factor type: BFlinearModel, JZS

# BF for CT quantity
lmBF_ctcmax / lmBF_covonly

Bayes factor analysis
--------------
[1] ctc_hour_max + cesd_t4 + NEG : 2.803794 ±0%

Against denominator:
  itsea_symptoms ~ NEG + cesd_t4 
---
Bayes factor type: BFlinearModel, JZS

# BF for infant vocalizations
mBF_cvcmax / lmBF_covonly

Bayes factor analysis
--------------
[1] cvc_hour_max + cesd_t4 + NEG : 4.723113 ±0%

Against denominator:
  itsea_symptoms ~ NEG + cesd_t4 
---
Bayes factor type: BFlinearModel, JZS

Create residualized symptoms variable for covariates in best fitting model

best_cov <- lm(
  itsea_symptoms ~
    scale(cesd_t4, scale = FALSE) +
    scale(NEG, scale = FALSE),
  data = d0
)

d0 <-
  d0 %>% 
  add_residuals(best_cov, var = "resid_symptoms")

Plot associations

AW consistency

d0 %>% 
  ggplot(aes(awc_prop, resid_symptoms)) +
  geom_point(size = 4, alpha = 1/2) +
  geom_smooth(method = "lm", size = 3, color = "black") +
  scale_x_continuous(
    breaks = seq.int(0, 1, .10),
    labels = scales::percent_format(accuracy = 1)
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 17),
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.position = "none"
  ) +
  labs(
    x = "Consistency of adult words in daily life",
    y = "Toddler symptoms of psychopathology\n(residuals)"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/aw_prop_symptoms.png",
  width = 7,
  height = 6
)

CT consistency

d0 %>% 
  ggplot(aes(ctc_prop, resid_symptoms)) +
  geom_point(size = 4, alpha = 1/2) +
  geom_smooth(method = "lm", size = 3, color = "black") +
  scale_x_continuous(
    breaks = seq.int(0, 1, .10),
    labels = scales::percent_format(accuracy = 1)
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 17),
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.position = "none"
  ) +
  labs(
    x = "Consistency of conversational turns in daily life",
    y = "Toddler symptoms of psychopathology\n(residuals)"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/ctc_prop_symptoms.png",
  width = 7,
  height = 6
)

Consistency

d0 %>% 
  dplyr::select(ctc_prop, awc_prop, resid_symptoms) %>% 
  rename(
    `CT consistency` = ctc_prop,
    `AW consistency` = awc_prop
  ) %>% 
  gather(measure, value, -resid_symptoms) %>% 
  ggplot(aes(value, resid_symptoms)) +
  geom_point(size = 4, alpha = 1/2) +
  geom_smooth(method = "lm", size = 3, color = "black") +
  scale_x_continuous(
    breaks = seq.int(0, 1, .15),
    labels = scales::percent_format(accuracy = 1)
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 17),
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.position = "none"
  ) +
  labs(
    x = "Consistency of language input in daily life",
    y = "Toddler symptoms of psychopathology\n(residuals)"
  ) +
  facet_grid(.~measure, scales = "free")

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/ctc_awc_prop_symptoms.eps",
  width = 10,
  height = 6,
  dpi = 1000
)

Infant vocalizations

d0 %>% 
  ggplot(aes(cvc_hour_max, resid_symptoms)) +
  geom_point(size = 4, alpha = 1/2) +
  geom_smooth(method = "lm", size = 3, color = "black") +
  scale_x_continuous(
    breaks = seq.int(0, 1000, 100)
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 17),
    panel.grid = element_blank(),
    plot.title = element_text(size = 18, hjust = .5),
    axis.title = element_text(size = 18),
    axis.text = element_text(size = 16),
    legend.position = "none"
  ) +
  labs(
    x = "Quantity of infant vocalizations in daily life",
    y = "Toddler symptoms of psychopathology\n(residuals)"
  )

ggsave(
  "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cv_max_symptoms.png",
  width = 7,
  height = 6
)

Elastic net regression with LENA variables and covariates

Prepare data

z_score <- function(x) {
  diff_mu <- x - mean(x, na.rm = T)
  sd <- sd(x, na.rm = T)
  diff_mu / sd
}

d0_en <-
  d0 %>% 
  dplyr::select(
    itsea_symptoms,
    ctc_prop,
    awc_prop,
    ctc_hour_max,
    awc_hour_max,
    cvc_hour_max,
    cesd_t4,
    NEG
  ) %>% 
  mutate_at(
    vars(itsea_symptoms:NEG),
    funs(z_score)
  ) %>% 
  na.omit()

predictors <- 
  d0_en %>% 
  dplyr::select(-itsea_symptoms) %>% 
  as.matrix()
  
symptoms <- z_score(d0_en$itsea_symptoms)

Run model over 100 increasingly relaxed lambda values

fit_net <- glmnet(predictors, symptoms, family = "gaussian")
plot(fit_net, label = TRUE)

Run leave-one-out cross validation to identify optimal lambda value

If nfolds is set to the sample size = leave-one-out CV https://sciphy-stats.com/post/2019-01-25-finalizing-glmnet-models/ Running the CV 100 times to minimize simulation error.

# run cross validation
set.seed(123)

lambdaMins <- c()

for (i in 1:100) {
  fit_cv <- cv.glmnet(
    predictors, 
    symptoms, 
    type.measure = "mse", 
    family = "gaussian", 
    alpha = .5,
    nfolds = 96,
    grouped = FALSE
  )
  lambdaMins <- cbind(lambdaMins, fit_cv$lambda.min)
}

# save lambda corresponding to minimal mse to a tibble
lambdaMins <- as_tibble(lambdaMins)

The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

lambdaMins <-
  lambdaMins %>%
  gather(run, lambdaMin, V1:V100) 

# calculate mean lambda value 
lambdaFinal <- mean(lambdaMins$lambdaMin)

Print coefficients for model with optimal lambda

fit_net_final <- glmnet(predictors, symptoms, family = "gaussian", alpha = .5, lambda = lambdaFinal)
coefs <- coef(fit_net_final) 

en_coef_tbl <- 
  tidy(coef(fit_net_final))  %>% 
  dplyr::select(-column) %>% 
  rename(
    Variable = row,
    `Elastic net estimate` = value
  ) %>% 
  mutate(
    Variable = recode(
      Variable,
      "(Intercept)" = "(Intercept)",
      "ctc_prop" = "CT consistency",
      "awc_prop" = "AW consistency",
      "ctc_hour_max" = "CT quantity",
      "cvc_hour_max" = "Infant vocalizations",
      "cesd_t4" = "Maternal depressive symptoms",
      "NEG" = "Infant negative affectivity"
    )
  ) %>% 
  mutate_at(vars(`Elastic net estimate`) , round, 2)

'tidy.dgCMatrix' is deprecated.
See help("Deprecated")'tidy.dgTMatrix' is deprecated.
See help("Deprecated")

Calculate EN model performance

# save model matrix
X <- model.matrix(itsea_symptoms ~ ., d0_en)

# get predicted values
y_hat <- matrix(X, ncol = 8) %*% as.matrix(coefs)

# calculate R-squared
r_sq <- 1 - (sum((y_hat[, 1] - symptoms)^2) / sum((symptoms - mean(symptoms))^2))

Evaluate model performance

for (i in 1000) {
  fit_rr_cv <- cv.glmnet(predictors, symptoms, type.measure = "mse", alpha = .5, family = "gaussian")
  lambdaMins <- cbind(lambdaMins, fit_rr_cv$lambda.min)
}

Additional analyses

Covarying maternal sensitivity

AW

sym_aw_ms <- 
  lm(
    scale(itsea_symptoms) ~ 
      scale(awc_prop) + 
      scale(cesd_t4) +
      scale(NEG) +
      scale(sens_fp) +
      scale(sens_r_sfp),
    data = d0 
  )


calc.relimp(sym_aw_ms)

Response variable: scale(itsea_symptoms) 
Total response variance: 1.019533 
Analysis based on 89 observations 

5 Regressors: 
scale(awc_prop) scale(cesd_t4) scale(NEG) scale(sens_fp) scale(sens_r_sfp) 
Proportion of variance explained by model: 24.51%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                          lmg
scale(awc_prop)   0.052732752
scale(cesd_t4)    0.028702483
scale(NEG)        0.134788698
scale(sens_fp)    0.021700415
scale(sens_r_sfp) 0.007145296

Average coefficients for different model sizes: 

                           1X         2Xs         3Xs         4Xs        5Xs
scale(awc_prop)   -0.24321526 -0.23740095 -0.23135613 -0.22535320 -0.2196113
scale(cesd_t4)     0.17364596  0.16915053  0.16557711  0.16284424  0.1608910
scale(NEG)         0.40773727  0.39466127  0.38145121  0.36835829  0.3556209
scale(sens_fp)    -0.19538884 -0.17194599 -0.14911310 -0.12717659 -0.1064187
scale(sens_r_sfp)  0.08106572  0.08337868  0.08447331  0.08445072  0.0834660

model_parameters(sym_aw_ms)

Parameter   | Coefficient |   SE |         95% CI |     t | df |      p
-----------------------------------------------------------------------
(Intercept) |       -0.04 | 0.10 | [-0.23,  0.15] | -0.44 | 83 | 0.662 
awc_prop    |       -0.22 | 0.10 | [-0.41, -0.03] | -2.26 | 83 | 0.026 
cesd_t4     |        0.16 | 0.09 | [-0.03,  0.35] |  1.69 | 83 | 0.094 
NEG         |        0.36 | 0.10 | [ 0.16,  0.56] |  3.54 | 83 | < .001
sens_fp     |       -0.11 | 0.10 | [-0.31,  0.10] | -1.04 | 83 | 0.300 
sens_r_sfp  |        0.08 | 0.09 | [-0.11,  0.27] |  0.88 | 83 | 0.382 

sym_ct_ms <- 
  lm(
    scale(itsea_symptoms) ~ 
      scale(ctc_prop) +
      scale(cesd_t4) +
      scale(NEG) +
      scale(sens_fp) +
      scale(sens_r_sfp),
    data = d0 
  )

calc.relimp(sym_ct_ms)

Response variable: scale(itsea_symptoms) 
Total response variance: 1.019533 
Analysis based on 89 observations 

5 Regressors: 
scale(ctc_prop) scale(cesd_t4) scale(NEG) scale(sens_fp) scale(sens_r_sfp) 
Proportion of variance explained by model: 25.17%
Metrics are not normalized (rela=FALSE). 

Relative importance metrics: 

                          lmg
scale(ctc_prop)   0.055000215
scale(cesd_t4)    0.029029047
scale(NEG)        0.140378673
scale(sens_fp)    0.020233226
scale(sens_r_sfp) 0.007072563

Average coefficients for different model sizes: 

                           1X         2Xs         3Xs        4Xs         5Xs
scale(ctc_prop)   -0.23739588 -0.23580976 -0.23431447 -0.2331960 -0.23268930
scale(cesd_t4)     0.17364596  0.16932768  0.16620912  0.1643043  0.16364810
scale(NEG)         0.40773727  0.39784264  0.38836046  0.3795487  0.37165793
scale(sens_fp)    -0.19538884 -0.16943354 -0.14329989 -0.1171948 -0.09132713
scale(sens_r_sfp)  0.08106572  0.08349367  0.08442111  0.0838533  0.08183789

model_parameters(sym_ct_ms)

Parameter   | Coefficient |   SE |         95% CI |     t | df |      p
-----------------------------------------------------------------------
(Intercept) |       -0.05 | 0.10 | [-0.24,  0.14] | -0.56 | 83 | 0.579 
ctc_prop    |       -0.23 | 0.10 | [-0.42, -0.04] | -2.43 | 83 | 0.017 
cesd_t4     |        0.16 | 0.09 | [-0.02,  0.35] |  1.73 | 83 | 0.087 
NEG         |        0.37 | 0.10 | [ 0.17,  0.57] |  3.73 | 83 | < .001
sens_fp     |       -0.09 | 0.10 | [-0.29,  0.11] | -0.89 | 83 | 0.374 
sens_r_sfp  |        0.08 | 0.09 | [-0.11,  0.27] |  0.87 | 83 | 0.389 

cor.test(d0$sens_fp, d0$itsea_symptoms)

    Pearson's product-moment correlation

data:  d0$sens_fp and d0$itsea_symptoms
t = -1.8456, df = 94, p-value = 0.0681
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.37348462  0.01401018
sample estimates:
      cor 
-0.187001 

cor.test(d0$sens_r_sfp, d0$itsea_symptoms)

    Pearson's product-moment correlation

data:  d0$sens_r_sfp and d0$itsea_symptoms
t = 0.79113, df = 92, p-value = 0.4309
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1224545  0.2801529
sample estimates:
       cor 
0.08220199 

Who is talking?

d0_follow <-
  d0 %>% 
  filter(
    lena_caregivers___1 != 0 |
      lena_caregivers___2 != 0 |
      lena_caregivers___3 != 0 |
      lena_caregivers___4 != 0 |
      lena_caregivers___5 != 0 |
      lena_caregivers___6 != 0 |
      lena_caregivers___7 != 0 |
      lena_caregivers___8 != 0 |
      lena_caregivers___9 != 0 |
      lena_caregivers___10 != 0
  ) %>% 
  rename(
    mom_present = lena_caregivers___1,
    dad_present = lena_caregivers___2,
    daycare_present = lena_caregivers___3,
    nanny_present = lena_caregivers___4,
    sibling_present = lena_caregivers___5,
    cousin_present = lena_caregivers___6,
    grandparent_present = lena_caregivers___7,
    friends_present = lena_caregivers___8,
    relative_present = lena_caregivers___9, #any other relative
    other_present = lena_caregivers___10
  )

d0_follow %>% 
  count(!is.na(percent_mother))

d0_follow %>%
  summarise_at(
    vars(percent_mother),
    funs(mean, median, sd, min, max), na.rm = TRUE
  )

cor.test(d0_follow$percent_mother, d0_follow$ctc_prop, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$ctc_prop
S = 81475, p-value = 0.0774
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.2065836 

cor.test(d0_follow$percent_mother, d0_follow$awc_prop, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$awc_prop
S = 72960, p-value = 0.4954
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
-0.08049238 

cor.test(d0_follow$percent_mother, d0_follow$ctc_hour_max, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$ctc_hour_max
S = 85137, p-value = 0.0248
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.2608283 

cor.test(d0_follow$percent_mother, d0_follow$awc_hour_max, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$awc_hour_max
S = 76045, p-value = 0.2841
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.1261794 

cor.test(d0_follow$percent_mother, d0_follow$cvc_hour_max, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$cvc_hour_max
S = 78635, p-value = 0.1613
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
-0.164533 

cor.test(d0_follow$percent_mother, d0_follow$itsea_symptoms, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$percent_mother and d0_follow$itsea_symptoms
S = 61876, p-value = 0.4786
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.0836531 

d0_follow %>%
  count(mom_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(dad_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(mom_present, dad_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(daycare_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(nanny_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(grandparent_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(friends_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(relative_present) %>% 
  mutate(
    per = n / sum(n)
  )

d0_follow %>%
  count(
    mom_present, 
    dad_present, 
    daycare_present, 
    nanny_present, 
    grandparent_present
  ) %>% 
  mutate(
    per = n / sum(n)
  ) %>% 
  arrange(n)

NA

d0_follow %>% 
  count(percent_mother) %>% 
  arrange(desc(n))

d0_follow <-
  d0_follow %>% 
  mutate(
    n_caregivers = 
      mom_present +
      dad_present +
      daycare_present +
      grandparent_present +
      nanny_present +
      friends_present +
      relative_present
  )

cor.test(d0_follow$ctc_prop, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$ctc_prop and d0_follow$n_caregivers
S = 58928, p-value = 0.09238
alternative hypothesis: true rho is not equal to 0
sample estimates:
     rho 
0.194426 

cor.test(d0_follow$awc_prop, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$awc_prop and d0_follow$n_caregivers
S = 61071, p-value = 0.154
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1651252 

cor.test(d0_follow$ctc_hour_max, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$ctc_hour_max and d0_follow$n_caregivers
S = 69248, p-value = 0.6472
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
0.05334703 

cor.test(d0_follow$awc_hour_max, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$awc_hour_max and d0_follow$n_caregivers
S = 46296, p-value = 0.001106
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.3671023 

cor.test(d0_follow$cvc_hour_max, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$cvc_hour_max and d0_follow$n_caregivers
S = 72778, p-value = 0.9653
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
0.005081446 

cor.test(d0_follow$itsea_symptoms, d0_follow$n_caregivers, method = "spearman")

Cannot compute exact p-value with ties

    Spearman's rank correlation rho

data:  d0_follow$itsea_symptoms and d0_follow$n_caregivers
S = 80673, p-value = 0.3767
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.1028426 

LENA recording correlations: duration, etc.

d0 %>% 
  dplyr::select(
    ctc_prop,
    awc_prop,
    ctc_hour_max,
    awc_hour_max,
    duration_analyzed,
    duration_total,
    ct_first_time,
    ct_last_time
  ) %>% 
  correlate(method = "spearman") %>% 
  fashion()

Correlation method: 'spearman'
Missing treated using: 'pairwise.complete.obs'

cor.test(d0$ctc_prop, d0$duration_analyzed, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_prop and d0$duration_analyzed
t = -0.8847, df = 98, p-value = 0.3785
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2805269  0.1093159
sample estimates:
        cor 
-0.08901346 

cor.test(d0$awc_prop, d0$duration_analyzed, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_prop and d0$duration_analyzed
t = -2.4278, df = 98, p-value = 0.01702
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.41518015 -0.04381767
sample estimates:
       cor 
-0.2381858 

cor.test(d0$ctc_hour_max, d0$duration_analyzed, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_hour_max and d0$duration_analyzed
t = 1.509, df = 98, p-value = 0.1345
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.04712287  0.33712968
sample estimates:
    cor 
0.15069 

cor.test(d0$awc_hour_max, d0$duration_analyzed, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_hour_max and d0$duration_analyzed
t = 1.3002, df = 98, p-value = 0.1966
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.06793023  0.31849685
sample estimates:
      cor 
0.1302255 

cor.test(d0$ctc_prop, d0$ct_first_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_prop and d0$ct_first_time
t = -1.7196, df = 98, p-value = 0.08865
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.35561061  0.02615038
sample estimates:
       cor 
-0.1711468 

cor.test(d0$awc_prop, d0$ct_first_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_prop and d0$ct_first_time
t = -0.20145, df = 98, p-value = 0.8408
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2159007  0.1767792
sample estimates:
        cor 
-0.02034532 

cor.test(d0$ctc_hour_max, d0$ct_first_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_hour_max and d0$ct_first_time
t = -2.1017, df = 98, p-value = 0.03814
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.38825834 -0.01173866
sample estimates:
       cor 
-0.2076779 

cor.test(d0$awc_hour_max, d0$ct_first_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_hour_max and d0$ct_first_time
t = -2.0029, df = 98, p-value = 0.04795
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.37992144 -0.00195945
sample estimates:
       cor 
-0.1983012 

cor.test(d0$ctc_prop, d0$ct_last_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_prop and d0$ct_last_time
t = -2.6758, df = 98, p-value = 0.008739
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.43505086 -0.06799646
sample estimates:
       cor 
-0.2609297 

cor.test(d0$awc_prop, d0$ct_last_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_prop and d0$ct_last_time
t = -2.4694, df = 98, p-value = 0.01526
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.41854994 -0.04788765
sample estimates:
       cor 
-0.2420292 

cor.test(d0$ctc_hour_max, d0$ct_last_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ctc_hour_max and d0$ct_last_time
t = -0.92389, df = 98, p-value = 0.3578
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2841549  0.1054191
sample estimates:
        cor 
-0.09292313 

cor.test(d0$awc_hour_max, d0$ct_last_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$awc_hour_max and d0$ct_last_time
t = -1.3196, df = 98, p-value = 0.19
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.32023976  0.06599806
sample estimates:
       cor 
-0.1321329 

cor.test(d0$ct_last_time, d0$ct_first_time, method = "pearson")

    Pearson's product-moment correlation

data:  d0$ct_last_time and d0$ct_first_time
t = 6.7365, df = 98, p-value = 0.000000001118
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4116546 0.6834769
sample estimates:
      cor 
0.5625842 

d0 %>% 
  ggplot(aes(duration_analyzed, awc_prop)) +
  geom_point() +
  geom_smooth(method = "lm")

d0 %>% 
  ggplot(aes(ct_first_time, awc_hour_max)) +
  geom_point() +
  geom_smooth(method = "lm")

d0 %>% 
  ggplot(aes(ct_first_time, ctc_hour_max)) +
  geom_point() +
  geom_smooth(method = "lm")

d0 %>% 
  ggplot(aes(ct_last_time, awc_prop)) +
  geom_point() +
  geom_smooth(method = "lm")

d0 %>% 
  ggplot(aes(ct_last_time, ctc_prop)) +
  geom_point() +
  geom_smooth(method = "lm")

CT and AW rate variables and symptoms

cor.test(d0$ctc_rate, d0$itsea_symptoms)

    Pearson's product-moment correlation

data:  d0$ctc_rate and d0$itsea_symptoms
t = -1.8434, df = 98, p-value = 0.06829
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.36631497  0.01384722
sample estimates:
       cor 
-0.1830688 

cor.test(d0$awc_rate, d0$itsea_symptoms)

    Pearson's product-moment correlation

data:  d0$awc_rate and d0$itsea_symptoms
t = -1.5658, df = 98, p-value = 0.1206
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.34214889  0.04146029
sample estimates:
       cor 
-0.1562301 

#write_csv(d0, "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/king_analyzed_data_20200528.csv")